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METHODS  FOR  CALCULATING  ATMOSPHERIC 
REFRACTION  AND  ITS  PERTURBATION 


INTRODUCTION 

In  obtaining  location  coordinates  on  the  Earth’s  surface  by  means  of  star 
observations,  errors  are  introduced  by  the  refraction  of  light  by  the  atmosphere 
because  the  actual  positions  of  the  stars  or  other  bodies  used  differ  from  their 
observed  positions.  A  knowledge  of  the  refractive  index  of  the  atmosphere  as  a 
function  of  position  permits  the  calculation  of  the  actual  ray  path  through  the 
atmosphere  by  numerical  integration  of  a  set  of  differential  equations  that  may  be 
used  to  extend  the  ray  step  by  step  through  the  atmosphere  to  obtain  the  true 
direction  of  the  ray  upon  emergence. 


In  this  report,  the  equations  and  numerical  procedures  for  performing  such 
a  ray  trace  are  derived  for  general  refractive  index  functions  including  lateral  re¬ 
fraction  effects,  and  explicit  forms  are  given  for  computation  of  the  refraction  error. 
Since  the  greater  variations  of  refraction  occur  with  altitude  and  because  the  de¬ 
tailed  state  of  the  atmosphere  is  usually  not  known  at  the  time  a  given  set  of  obser¬ 
vations  is  made,  it  is  desirable  (and  sufficient  for  many  purposes)  to  consider  the 
refractive  index  to  vary  only  with  the  radial  coordinate  and  to  neglect  the  effects 
due  to  lateral  variations.  For  this  spherically  symmetrical  case,  the  refraction  error 
is  present  only  for  the  altitude  angle  and  is  reduced  to  a  quadrature  over  the  radial 
distance  or  height.  The  strongest  part  of  the  error  can  be  considered  to  originate 
from  average  atmospheric  profiles  that  can  be  approximated  by  simple  formulas 
and  several  of  these  standard  atmospheric  representations  that  have  commonly  been 
used  are  presented.  The  perturbation  of  the  refraction  error  due  to  variation  of  the 
atmosphere  from  the  standard  atmosphere  leads  to  a  further  correction  that  is 
derived  and  given  in  terms  of  another  quadrature  formula. 
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RAY  TRACING  IN  SPHERICAL  COORDINATES 

According  to  Fermat’s  principle  (also  referred  to  as  the  principle  of  least 
time),1  the  ray  joining  any  two  arbitrary  points,  Pj  and  P2 ,  is  determined  by  the 
condition  that  its  optical  length 


(1) 


be  stationary  as  compared  with  the  optical  lengths  of  arbitrary  neighboring  curves 
joining  Pj  and  P2 .  If  the  refractive  index  n  is  considered  to  be  a  given  smooth 
continuous  function  of  position  and  the  location  along  the  path  is  given  in  terms  of 
a  parameter  t,  then  an  actual  ray  path  must  furnish  an  extremum 


p 

5  f  2  n(r,0,0)S(r,0,;,0,*)dt  =  0  (2) 


where  spherical  coordinates  are  indicated  with 

=  S(r,  0,r,  0,0)  =  ^~t262~+~t2~ sin2^2""  (3) 


and  where  the  dots  indicate  differentiation  with  respect  to 't.  The  partial  derivatives 


as 

r(0*2  +  sin2  0  02) 

as 

• 

r 

dr 

S 

”  a? 

~s~ 

3S 

r2  sin  0  cos  0  <p 2 

as 

r20 

30 

S 

~a£ 

S 

as 

b<t> 

0 

as 

’  w 

r2  sin2  0  0 
S 

*H.A.  Buchdahl,  An  Introduction  to  Hamiltonian  Optics,  Cambridge  University  Press,  1970. 
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the  differential  system  for  the  rays  is  simplified  by  eliminating  the  radicals  appearing 
in  S  above. 


(9) 


-  (nr)  -  — - nr(0  +  sin  6  <f>2)  =  0 

ds  dr 


—  (nr2  0)  — 777 —  nr2  sin  0  cos  6  <p2  =  0 
ds  ad 


d  ,  .  ,  •  dn 

—  ( nr 2  sin2  0  0)  -  -r—  =  0 
ds  d0 


If  a  canonical  system  of  variables  is  introduced  where 


Pr  =  nr 
Pfl  =  nr20 
p0  =  nr2  sin2  0  (p 


(10) 


the  corresponding  first-order  differential  system  is  easily  put  in  normal  form. 

2 


Pr 


1  ,  p*  dn 

=  - (P.2  +  )  +  — 


nrJ 


sur  0 


) 


dr 


.  cos  0  p*2  dn 

p.  =  - + 


nr2  sin3  0  d0 


PA  * 


dn 

dtp 


(11) 


r  = 


Pt_ 

n 


0  = 


nr 


<p  = 


nr2  sin2  0 


This  system  is  suitable  for  numerical  integration  by  many  standard  methods  in¬ 
cluding  the  Runge-Kutta  method.  The  equations  are  not  completely  independent 
but  are  inter-related  by  the  implicit  relationship  from  equations  3  and  8 

r2  +  r2  02  +  r2  sin2  9  Q2  =  1  (12) 

which  requires  that  the  sum  of  the  squares  of  the  local  direction  cosines  of  the 
tangent  to  the  ray  at  any  point  be  unity.  This  enables  the  integration  to  be  initiated 
from  a  knowledge  of  position  coordinates  and  two  angles  sighted  along  a  ray;  for 
example,  altitude  and  azimuth  angles.  It  can  also  facilitate  the  change  of  inde¬ 
pendent  variable  from  s  to  one  of  the  coordinates  if  desired;  for  example,  if  it 
is  desired  to  increment  the  radial  distance  r  in  fixed  predetermined  amounts.  In 
such  a  case,  the  six  equations  given  by  11  are  reduced  to  five.  For  a  general  inte¬ 
gration  with  s  as  the  independent  variable,  the  initial  conditions  consist  of  the 
coordinates  r,  0,  0  and  the  direction  cosines  ar,  ag ,  a0,  related  to  the  con¬ 
jugate  variables,  as  follows. 


Pr  =  n«r 


a°  r  ds 


P«  =  nrQ„ 


•  a 

a.  =  r  sinfl  — — 
*  ds 


P*  =  nr  sin  da^ 


The  altitude  angle  a  and  azimuth  angle  A  are  given  by 


where  the  ambiguity  of  sign  must  be  rectified  to  conform  with  the  spherical  co¬ 
ordinates,  since  various  defining  conventions  are  used  for  azimuth.  By  using  the 
identity 


K: 


<*r2  +  “92  +  %2  =  (15) 

it  is  easy  to  obtain  the  direction  cosines  in  terms  of  altitude  and  azimuth. 

0ir  =  sin  a 

ol$  =  cos  a  cos  A  (16) 

a .  =  +  cos  a  sin  A 

<p  — 


S 

ATMOSPHERIC  REFRACTION  INCLUDING 
LATERAL  REFRACTION 


* 

■■ 


dn  dm  dn 

Assuming  the  quantities  dr.  qq  are  known  functions  of  position, 

one  may  trace  a  ray  up  through  the  atmosphere  by  using  equations  1 1  for  any 
starting  location  r0 ,  dQ ,  4>0  and  direction  aIo ,  ag  q  .  Assuming  the  initial 
altitude  angle  is  great  enough  that  atmospheric  ducting  ant?  subsequent  return  of 
the  ray  does  not  occur,  the  ray  eventually  will  emerge  from  the  atmosphere  at 
some  location  rf,  0f,  <t>{  with  local  direction  coordinates  arf>  a6{,  a0f.  In 
order  to  determine  the  amount  of  bending  of  the  ray,  the  transformation  of  the 
final  direction  coordinates  back  into  the  initial  frame  must  be  known.  This  trans¬ 
formation  will  now  be  obtained.  For  a  general  position  vector  IT  given  in  rec¬ 
tangular  components  but  expressed  in  spherical  coordinates 

-*•  A  A  A 

R  =  i  r  sin  d  cos  <p  +  j  r  sin  6  sin  (j>  +  k  r  cos  9.  (17) 


U 
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A  A  A 

A  local  reference  frame  of  unit  vectors  r,  0 ,  0  may  be  defined  by 


3R  /  I  _3R 
dr  /  |  3r 


A  A  A 

=  i  sin  0  cos  0  +  j  sin  0  sin  0  -  k  cos  0 


3R 

30 


/ 


3  R  A  A  A 

-rr^-  =  i  cos  0  cos  0  +  j  cos  0  sin  0  +  k  sin  0 

ou  (18) 


0A=  J*  / 

9  30  / 


3R 

30 


A  A 

=  -  i  sin  0  +  j  cos  0 


If  the  unit  direction  vector  of  the  ray 


a*  =  a  i  +  a  0  +  cr  0 

r  8  0  Y 


(19) 


is  expressed  in  terms  of  the  initial  frame,  the  components  are  found  to  depend 
on  the  cosines  of  angles  between  the  initial  and  current  frame  vectors 

t(r  .  ?0)at  +  (0  •  P0)«8  +  (0  .  r0)a#l 

+  •  *„>«,  +  (i  ' *  »oK]  (20) 

#0  $0)ar  +  (0  •  0o)a8  +  (0  •  0O)%  J 

where  the  direction  cosines  involved  are  readily  obtained  from  equations  1 8  applied 
at  the  initial  and  current  positions  (The  prime  added  to  a'  is  to  avoid  confusion 
with  the  starting  direction  ocQ ). 


+ 


+ 
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r  •  =  sin  0  sin  0O  cos  (0  -  0O)  +  cos  0  cos  0O 

r  .$0  =  sin  0  cos0o  cos(0-0o)  -  cos0  sin  0O 

r  •  $0  =  sin0sin(0-0o) 

A  /\ 

0  •  r0  =  cos0  sin0o  cos(0-0o)  -  sin0  cos0o  (21) 

$  •  =  cos  0  cos  0O  cos  (0  -  0O  )  +  sin  0  sin  0O 

$  •  0  =  cos  0  sin  (0  -  0O  ) 

A  A 

0  *ro  =  -  sin  0O  sin  (0  -  0O) 

A  A 

0  4  0o  =  -cos0Q  sin  (0  —  0O) 

$  *$o  =  cos  (0  -0O) 


By  applying  equations  20  and  21  to  the  emerging  ray  direction  af,  the  com¬ 
ponents  referred  to  the  initial  frame  can  be  expressed  in  matrix  form  as  given  by 
equation  22. 


“,'r 

sinflf  sinflo  cos(0r-0o)  +  cos0(cos0o 

cos  0,  sin  0^  cos(0r  -  0o)  -  sin  0r  cos  0o 

-sin0osin(0f-*o7 

“if 

= 

sin  0f  cos  0o  cos(0f  -  4a)  -  cos  0f  sin  8g 

cos  0f  cos  0g  cos(0r  -  4g )  +  sin  0f  sin  8g 

-cos  0o  sin(#f  -  ) 

sin  0f  sin  (4,  -40) 

cos  0f  sin  (4,  -  4Q ) 

cos  («,  -*„> 

"♦< 

It  was  found  that  the  transformation  matrix  could  be  factored  as 


“r'f 

sin  0 

O 

cos  0 

O 

0 

cos  0 

To 

sin  0 

O 

0 

“of 

= 

COS  dQ 

-sin  0O 

0 

0 

0 

1 

0 

0 

1 

-sin  0O 

COS  0O 

0 

cos  0f  sin  0f  0 

sin  <j>f  -cos  <p{  0 

0  0  1  ] 
(23) 


sin0f  cos0f  0 

“rf 

0  0-1 

cos  0f  -sin  0f  0 

or  alternatively  in  the  form  given  by  equation  24  as  probably  the  most  convenient 
for  computation. 


“r'f 

sin  0O  cos  0O  0 

cos (0f-0o)  sin  (0f  -  0O )  0 

= 

cos  0O  -sin  0O  0 

0  0  1 

%( 

0  0  1 

si n(0f-0o)  -cos(0f-0o)  0 

From  equation  14,  the  vertical  refraction  correction  is  given  by 

a/  -  a*  =  arccos  a7  -  arccos a  (25) 

to  n  ro  v  7 

and  the  lateral  refraction  error  by 


a' 


a. 


A/ -A  =  arctan  (+  0f— )  -  arctan  (±  — ^ — ) 

10  § 


a 


8  f 


a 


*o 


a7,  a  -a.  a7 


=  +  arctan  (- 


0f  9  ( 


$0  ®  f 


a7  a.  _a  a7 

®f  ®0  00  0f 


where  the  sign  again  depends  on  the  convention  used  for  azimuth. 


(26) 


THE  SPHERICALLY  SYMMETRICAL  CASE 

For  the  case  where  the  refractive  index  depends  only  upon  r,  equation  9 
becomes 


iL<n?)  -  --  -  nr(02  +  sin2  0  02)  =  0 
as  or 


-^-(n  r20)  -  nr2  sin  0  cos  0  <p2  =  0  (27) 

ds 


nr2  sin2  0  <p  =  C. 


where  an  integral  has  been  found  for  the  last  equation.  The  coordinate  system  may 
be  chosen  so  that  initially  -g*-  =  0.  Then,  Cj  =  0  and  -*j|  vanishes  indentically 

<t>  =  <j>0  =  constant  (28) 


and  the  problem  is  reduced  to  two  dimensions.  From  the  fact  that  <p  =  0,  the 
second  equation  of  27  becomes  integrable: 


nr 


2 


dfl  = 
ds 


c2 


(29) 


dfl 

Inserting  the  resultant  value  for  ds  into  the  first  equation  of  27  (together  with 
0  =  0)  yields  the  following  relationship: 


d  dr  3n  C  2 

- — (n- — )-  - — - \- 

ds  ds  3r  nr-* 


(30) 


By  multiplying  by  n  and  using  the  relationship 


_d^ 

ds 


dr  d 
ds  dr, 


n  ^  (  dr  n  dn 

ds  dr  ds  dr 


and  integrating  yields 


(n-£-)2  -  n2 
ds 


7^  -  C3 


(31) 


(32) 


C22  dfl  7 

If  -j-  is  replaced  by  (nr  - V  from  equation  29,  it  is  found  that 

r  ds 


n2  [(— )2  +  (r— )2  -  1]  = 


ds 


ds 


'3‘ 


(33) 
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The  bracketed  quantity  must  vanish  because  arc  length  s  is  defined  by  ds  = 
\J  dr2  +  r2  d02  and  hence  C3  =  0.  It  then  follows  that  equations  29  and  30 
(in  r  and  6)  are  not  independent.  As  a  matter  of  convenience,  equation  29  will 
be  used  and  the  geometrical  relations  between  dr,  ds  and  dd  will  be  exploited. 


As  demonstrated  in  figure  1,  a  star  is  observed  at  the  apparent  position  A }  given 
by  angle  4>0. 


FIGURE  1 .  Geometrical  Parameters  for  the  Atmospheric  Ray  Path. 
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If  no  atmosphere  were  present,  the  actual  position  A2  would  coincide  with  . 
For  a  ray  traveling  in  the  reverse  direction  and  emanating  at  the  surface  at  angle 
4>0  in  a  medium  with  variable  refractive  index,  the  ray  path  is  curved  and  its  in¬ 
clination  V/  at  r  is  given  by 


nr(r  - )  =  nr  cos  ^  =  C  (34) 

ds 

from  equation  29  where  the  constant  is  determined  from  the  initial  values  of 
n,  r  and  ip. 


C  =  n0ro  cos*0  (35) 

After  passing  through  the  region  of  variable  index,  the  ray  will  emerge  at  rf,  0{ 
in  the  direction  ip{  toward  A2.  By  the  optical  Principle  of  Reversibility,  an 
object  at  A2  would  be  observed  to  have  elevation  \p0 ;  whereas,  if  the  refractive 
index  were  constant  (atmosphere  removed),  the  true  elevation  would  be  angle 
a .  As  layers  of  variable  refractive  index  are  added  in  the  reversed  ray  system,  a 
would  change  and  so  in  this  inbedded  sense  can  be  regarded  as  a  function  of  r. 


From  figure  1 , 


a  =  i//  -  0 
or 

da  _  dip  d  0 

dr  dr  dr 


(36) 

(37) 


and  a  may  be  determined  by  integrating  equation  37.  By  rewriting  equation 
34  in  the  form 


and  using  the  fact  that  sin  i ]/  =  - ,  an  expression  for  - 

ds  dr 


sin  t|/ 


dd 


dr  nr2 


By  differentiating  equation  34  in  the  form 


cos  \J/  = 


nr 


d  d/ 

an  expression  containing  -  is  obtained. 

dr 


sin  i p 


d\p 


C  C  dn 
+ 


dr  nr^  n2r  dr 


Combining  equations  39  and  41 


sin  \p 


d\p  dd 


dr  dr 


C  dn 
n2r  dr 


and  using  the  fact  that 


sin 


=  \J  1  -  COS2  \l>  =  y/ 1  - 


an  expression  for  — —  is  readily  found. 


If  this  expression  is  integrated  by  parts  from  r0  to  rf)  the  value  for  the  obser¬ 
vational  error  is  found: 


r  dn 


(45) 


6a 


where  aQ  =  \p0,  iQ,  nQ  and  af ,  r{,  nf  are  initial  and  final  values,  and  C  is 
given  by  equation  35. 


r  (47) 

1 


For  determining  5a  by  numerical  integration,  equation  45  should  be  preferable 
to  equation  47  by  virtue  of  its  simplicity  and  certainly  a  need  to  carry  fewer 

significant  figures.  It  can  be  easily  evaluated  with  the  trapezoidal  rule,  using  a 

dn 

linear  interpolation  for  — —  For  higher  degree  approximations,  standard  spline 

dr 

methods  are  suggested.  Although  it  can  be  integrated  by  quadrature  formulae 
(e.g.  Newton-Cotes),  equation  47  appears  to  offer  no  distinct  advantage. 
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NUMERICAL  INTEGRATION 


The  equation  for  numerical  integration  with  the  trapezoidal  rule  using  linear 
cerpolation  for  the  derivative  can  be  copied  directly  from  equation  45  by  in¬ 
spection.  If  the  interval  [r0,  rf]  is  broken  into  N  sub-intervals  such  that 

r0  <...<  rj  <  rj+i  <...<  rN  =  rf  (48) 


and  where 


n,  =  nOj) 


(49) 


"i  =  ni  "  ni-i 

the  value  of  the  integral  is  given  approximately  by 


6a  =  Y, 


f.  +  f. 

j  j-i 


"j 


(50) 


where  C  =  n0r0  cos  <|/0  as  previously  defined.  This  form  enables  one  to  use 
unequal  intervals  that  may  be  necessary  for  some  sets  of  data. 


ATMOSPHERIC  MODELS 


In  order  to  perform  an  integration  with  equations  50  or  11,  one  must 
know  the  refractive  index  as  a  function  of  position.  This  is  a  function  of  temper¬ 
ature,  pressure,  and  composition.  Older  references  such  as  Newcomb  give  the  re¬ 
fractive  index  in  terms  of  density  p  by2 

n2  -  1  =  2Cp  (51) 

where  the  value  C  =  0.22607  gives  good  agreement  with  the  well-known  Poulkova 
tables  (excerpted  by  Tricker3).  For  n  close  to  1,  the  expression  n2  -  1  may 
be  replaced  by  2(n  -  1 ). 

n  =  1  +  Cp  (52) 

This  expression  is  said  to  give  still  better  agreement  with  measured  refraction.4 


The  actual  values  of  temperature,  pressure,  and  composition  of  the  atmos¬ 
phere  vary,  depending  on  position  and  time,  and  follow  no  simply  expressible 
physical  law  exactly.  Some  tables  giving  detailed  values  of  temperature  and  pressure 
as  a  function  of  altitude  are  available.5  Hotine  cites  the  following  formula  that  was 
adopted  by  the  International  Association  for  Geodesy  in  I960:6 


Mg  ~  1  /  p  \  55  X  10'7  e 

1  +  at  \  760  /  1  +  at 


(53) 


y 

Simon  Newcomb,  Compendium  of  Spherical  Astronomy,  New  York:  Dover  Publications.  1960. 
3R.A.R.  Tricker,  Introduction  to  Meteorological  Optics,  Elsovier.  1971. 

4Simon  Newcomb,  Compendium  of  Spherical  Astronomy,  New  York:  Dover  Publications.  1960. 
5£/.S.  Standard  Atmosphere,  1962,  National  Aeronautics  and  Space  Administration. 

6Martin  Hotine,  Mathematical  Geodesy,  ESSA  Monograph  2.  1969. 
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where 


p  =  actual  refractive  index 

juG  =  group  refractive  index  calculated  from 

(PG  -  1)  X  107  =  2876.04  +  (3  X  16.288)X~2  +  (5  X  0.136)X~ 

t  =  temperature  of  air  in  °C 

p  =  total  atmospheric  pressure  in  mm.  Hg. 

e  =  partial  pressure  of  water  vapor  content  in  mm.  Hg.  run.  Hg. 

a  =  temperature  coefficient  of  refractivity  of  air  ur 

(or  the  coefficient  of  thermal  expansion), 

(0.003661) 

X  =  wavelength  in  microns 

An  equivalent  formula  is  given  by  Bomford.7 


Several  formulas  for  air  density  p  as  a  function  of  altitude  h  are  given  by 
Newcomb.8 


h 


P  =  Px  e'Trr 

(54) 

h  \ 

p  =  (l- T 

— ) 

(55) 

hi  / 

>kh 

p  -  Pt  e"BY 

(56) 

7 

G. Bomford,  Geodesy,  Oxford  University  Press.  1971. 

O 

Simon  Newcomb,  Compendium  of  Spherical  Astronomy,  New  York:  Dover  Publications.  1960. 


In  equation  63,  the  full  variation  is  obtained  for  e  =  1,  and  the  conditions  that 
the  first-order  term  give  a  good  representation  of  the  corresponding  variation  in 
J  are 


em  (r) 

n(r) 

(64) 

dm 

dn 

dr 

dr 

where  e  has  been  carried  in  equation 

64 

mainly  for  purposes  of  identification. 

By  differentiation  and  by  a  considerable  amount  of  algebraic  manipulation, 
the  following  identity  may  be  obtained: 


and  this  is  useful  in  further  simplifying  the  form  of  the  second  integral  of  equation 
63. 


mr  irf 


+  0(e2) 


(66) 


By  assuming  the  value  of  m  to  vanish  at  the  endpoints, 

m(r0)  =  m(rf)  =  0  (67) 

the  quantity  bracketed  in  equation  66  will  also  vanish.  For  the  upper  endpoint, 
this  is  a  reasonable  assumption  since  the  refractive  index  should  assume  the  value 
for  vacuum  and  variation  or  perturbation  is  not  reasonable.  For  the  lower  endpoint, 
it  is  necessary  on  practical  grounds,  since  any  variation  of  refractive  index  will 
disturb  the  value  of  C  (initial  condition)  used  throughout  the  entire  range  of 
integration. 


The  first  term  of  equation  66  is  the  unperturbed  error.  The  integral  in  the 
second  term  is  known  as  the  variational  or  functional  derivative  of  J  of  first  order. 
Taking  e  =  1  and  ignoring  higher  order  terms  yields  the  first  order  or  linear  per¬ 
turbation  of  J : 
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and  where  C0  is  given  by  equation  52  in  order  to  express  the  variation  of  re¬ 
fractive  index  in  terms  of  the  variation  in  density.  The  end  values  of  j  do  not 
appear  in  equation  70  because  the  variation  in  refractive  index  is  taken  to  be  zero. 
Explicit  retention  of  the  factors  Ar^  enables  one  to  use  unequal  intervals  that 
may  be  necessitated  by  available  data. 
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CONCLUSIONS 


By  using  Fermat’s  Principle  and  the  Euler  equations,  a  system  of  differential 
equations  for  ray  tracing  in  spherical  coordinates  was  derived  for  media  with  con¬ 
tinuously  variable  refractive  index.  This  system  was  reduced  to  first  order  in  a  form 
suitable  for  numerical  integration.  The  transformation  from  the  final  reference 
frame  to  the  initial  frame  was  derived  and  expressions  for  the  true  altitude  and 
azimuth  angles  were  derived.  The  spherically  symmetrical  case  was  discussed  and 
the  integral  for  the  error  in  altitude  angle  was  obtained  along  with  the  corresponding 
numerical  formulas.  Some  standard  atmospheric  profiles  and  relationships  between 
refractive  index  and  density  given  in  the  literature  were  described.  A  perturbation 
equation  was  derived  for  treating  the  error  in  altitude  angle  due  to  departures  of 
the  actual  refractive  index  profile  from  that  of  a  standard  reference  atmosphere. 
The  corresponding  numerical  integration  formulas  were  also  obtained. 
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